n-Body Simple Harmonic Oscillators

In a stark contrast to gravitational systems, n-body simple harmonic oscillators are always soluble. This is because the differential equations governing their motion are linear no matter how many interactions are included, so the solution is a matter of matrix algebra. It is straightforward to evaluate the exact solution for an arbitrary number of bodies in an arbitrary number of dimensions.

Consider an n-body system with equal masses and frequencies, since the solution is particularly elegant. The Lagrangian for pairwise interaction in d spatial dimensions is

$L = \frac{m}{2} \sum_{i = 1}^{n} \sum_{k = 1}^{d} ({\overset{\cdot}{x}}_{i}^{k})^{2} - \frac{m ω^{2}}{2} \sum_{i = 1}^{n - 1} \sum_{j = i + 1}^{n} \sum_{k = 1}^{d} (x_{i}^{k} - x_{j}^{k})^{2}$

where the upper indices label spatial dimensions as is done in general relativity. Evaluating the Lagrange equations

$\frac{d}{d t} \frac{\partial L}{\partial {\overset{\cdot}{q}}_{i}} = \frac{\partial L}{\partial q_{i}}$

for any one of the variables gives simply

${\overset{\cdot \cdot}{x}}_{i}^{k} = - ω^{2} \sum_{j \neq i} (x_{i}^{k} - x_{j}^{k})$

which are linear equations as stated at the outset. The overall sign on the right-hand side can be confirmed by expanding each individual square in the potential.

These equations are independent of the number of spatial dimensions in the sense that the same form holds for each Cartesian coordinate individually. Since a general solution for one coordinate holds for all, apart from initial conditions, the upper index can be dropped temporarily.

The right-hand side of these equations contains n−1 terms for the variable on the left-hand side and one term for each of the remaining variables. The set of equations for each Cartesian coordinate thus has a matrix on the right-hand side which is all ones apart from the diagonal:

$\frac{d^{2}}{d t^{2}} [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}] = ω^{2} [\begin{array}{c} 1 - n & 1 & \dots \\ 1 & 1 - n & \dots \\ ⋮ & ⋮ & ⋱ \end{array}] [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}] = ω^{2} M [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}]$

The general behavior of this matrix can be understood by splitting it into two matrices, one which is all ones and the other proportional to the identity matrix:

$M = J - n I J = [\begin{array}{c} 1 & 1 & \dots \\ 1 & 1 & \dots \\ ⋮ & ⋮ & ⋱ \end{array}]$

A matrix of all ones has the interesting property

$J^{2} = n J$

which is easy to understand because multiplying any column with any row simply gives the dimension of the matrix. From this one can determine

$M^{2} = (J - n I)^{2} = n J - 2 n J + n^{2} I = - n (J - n I) = - n M$

Multiplying repeatedly by M on both sides gives the general formula

$M^{k} = (- n)^{k - 1} M$

while taking the square root of both sides gives

$\sqrt{M} = \pm \frac{i}{\sqrt{n}} M$

consistent with the general formula. This last expression is needed for writing down the general solution to the system of equations for each Cartesian coordinate,

$[\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}] = exp (ω \sqrt{M} t) [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}]_{0} = exp (\pm \frac{i ω t}{\sqrt{n}} M) [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}]_{0}$

where the subscripts indicate initial values. Evaluating the exponential factor is straightforward:

$\begin{array}{l} exp (\pm \frac{i ω t}{\sqrt{n}} M) = I + \sum_{k = 1}^{\infty} (\pm \frac{i ω t}{\sqrt{n}})^{k} \frac{M^{k}}{k!} \\ = I + \sum_{k = 1}^{\infty} (\pm \frac{i ω t}{\sqrt{n}})^{k} \frac{(- n)^{k - 1}}{k!} M = I - \frac{M}{n} \sum_{k = 1}^{\infty} \frac{(\mp i \sqrt{n} ω t)^{k}}{k!} \\ = I - \frac{M}{n} [exp (\mp i \sqrt{n} ω t) - 1] = \frac{J}{n} - \frac{M}{n} exp (\mp i \sqrt{n} ω t) \end{array}$

In terms of real functions one then has

$\begin{array}{l} cos (\frac{ω t}{\sqrt{n}} M) = \frac{J}{n} - \frac{M}{n} cos (\sqrt{n} ω t) \\ sin (\frac{ω t}{\sqrt{n}} M) = \frac{M}{n} sin (\sqrt{n} ω t) \end{array}$

so that the general solution can be written

$\begin{array}{l} [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}] = cos (\frac{ω t}{\sqrt{n}} M) [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}]_{0} - \frac{1}{\sqrt{n} ω} sin (\frac{ω t}{\sqrt{n}} M) [\begin{array}{c} {\overset{\cdot}{x}}_{1} \\ {\overset{\cdot}{x}}_{2} \\ ⋮ \end{array}]_{0} \\ = [\frac{J}{n} - \frac{M}{n} cos (\sqrt{n} ω t)] [\begin{array}{c} x_{1} \\ x_{2} \\ ⋮ \end{array}]_{0} - \frac{M}{n^{3 / 2} ω} sin (\sqrt{n} ω t) [\begin{array}{c} {\overset{\cdot}{x}}_{1} \\ {\overset{\cdot}{x}}_{2} \\ ⋮ \end{array}]_{0} \end{array}$

Although the idempotent matrix M is not invertible, one can still multiply both sides of the initial condition for velocity to determine the proper relationship to input values. Note that the entire solution is divided by the number of bodies.

It should be stressed that this solution holds for all Cartesian coordinates and any number of interacting bodies. The system with equal masses and frequencies clearly has a tremendous amount of symmetry.

As a check of the general solution, consider two bodies interacting in one dimension, which is immediately separable:

$\begin{array}{l} {\overset{\cdot \cdot}{x}}_{1} = - ω^{2} (x_{1} - x_{2}) \\ {\overset{\cdot \cdot}{x}}_{2} = - ω^{2} (x_{2} - x_{1}) \end{array} \to \begin{array}{l} {\overset{\cdot \cdot}{x}}_{1} - {\overset{\cdot \cdot}{x}}_{2} = - 2 ω^{2} (x_{1} - x_{2}) \\ {\overset{\cdot \cdot}{x}}_{1} + {\overset{\cdot \cdot}{x}}_{2} = 0 \end{array}$

The factor here of two represents a reduced mass. The solutions to these equations are

$\begin{array}{l} x_{1} - x_{2} = (x_{1, 0} - x_{2, 0}) cos (\sqrt{2} ω t) + \frac{{\overset{\cdot}{x}}_{1, 0} - {\overset{\cdot}{x}}_{2, 0}}{\sqrt{2} ω} sin (\sqrt{2} ω t) \\ x_{1} + x_{2} = x_{1, 0} + x_{2, 0} + ({\overset{\cdot}{x}}_{1, 0} + {\overset{\cdot}{x}}_{2, 0}) t \end{array}$

The individual variables are then determined by adding and subtracting these two expressions. Comparing the result to the general solution indicates that it represents the same motion but with the center of mass unmoving, since there is no term linear in time. This check also indicates the source of the overall denominator of the number of bodies.

And now to see it in action! The frequency is taken as unity because it merely represents a rescaling of the temporal variable. Initial positions and velocities are randomized about the origin, with the velocities adjusted to keep the center of mass at the origin. In two dimensions the evolution of the system looks like this:

The size of the bodies has been increased for legibility, but one can imagine them as points that do not collide. The reset button starts the visualization again with random initial conditions.

Since the motion is completely soluble in terms of simple functions, it is clearly precisely periodic. There is no chaos possible in this system, no matter how many bodies are involved. Quite different from gravitation in this respect!

In three dimensions the evolution of the system looks like this:

Again, the bodies can be imagined as point spheres that do not collide. The coordinate axes are included to provide some perspective for the scene. Apart from the added spatial coordinate, this visualization is essentially the same as the last.

Although this system is simple enough to be exactly soluble, one would still expect it would have something to say about n-body problems for other pairwise potentials. To be continued...